Characterization of a prognostic model for lung squamous cell carcinoma based on eight stemness index-related genes

Background Cancer stem cells (CSCs) are implicated in cancer progression, chemoresistance, and poor prognosis; thus, they may be promising therapeutic targets. In this study, we aimed to investigate the prognostic application of differentially expressed CSC-related genes in lung squamous cell carcinoma (LUSC). Methods The mRNA stemness index (mRNAsi)-related differentially expressed genes (DEGs) in tumors were identified and further categorized by LASSO Cox regression analysis and 1,000-fold cross-validation, followed by the construction of a prognostic score model for risk stratification. The fractions of tumor-infiltrating immune cells and immune checkpoint genes were analyzed in different risk groups. Results We found 404 mRNAsi-related DEGs in LUSC, 77 of which were significantly associated with overall survival. An eight-gene prognostic signature (PPP1R27, TLX2, ANKLE1, TIGD3, AMH, KCNK3, FLRT3, and PPBP) was identified and used to construct a risk score model. The TCGA set was dichotomized into two risk groups that differed significantly (p = 0.00057) in terms of overall survival time (1, 3, 5-year AUC = 0.830, 0.749, and 0.749, respectively). The model performed well in two independent GEO datasets (p = 0.029, 0.033; 1-year AUC = 0747, 0.783; 3-year AUC = 0.746, 0.737; 5-year AUC = 0.706, 0.723). Low-risk patients had markedly increased numbers of CD8+ T cells and M1 macrophages and downregulated immune checkpoint genes compared to the corresponding values in high-risk patients (p < 0.05). Conclusion A stemness-related prognostic model based on eight prognostic genes in LUSC was developed and validated. The results of this study would have prognostic and therapeutic implications. Supplementary Information The online version contains supplementary material available at 10.1186/s12890-022-02011-0.

theory, CSCs are implicated in tumor initiation, growth, and metastasis [8,9]. CSCs are major contributors to the resistance to conventional therapies, tumor recurrence, and poor prognosis. Therefore, targeting CSCs offers a new approach to developing efficient therapies and improving outcomes [10,11]. An increasing evidence shows that CSCs have prognostic value and could serve as potential prognostic biomarkers in various cancers, including lung cancer [12][13][14]. Recently, Liao et al. identified key cancer stemness-related genes implicated in LUSC through integrated bioinformatics analysis [15]. Similarly, Qin et al. screened LUSC mRNA-related hub genes through bioinformatics and concluded that these genes may serve as therapeutic targets for inhibiting the LUSC stem cell properties [16]. Nevertheless, there is a lack of research on the application of cancer stemnessrelated genes as prognostic tools in LUSC.
Many studies have been conducted on the derivation of gene signatures as a way of determining prognostic potential; for example, Giannos et al. identified prognostic genetic biomarkers for cell lung cancer progression through comprehensive bioinformatics analysis [17], and Wu et al. identified hub genes and important KEGG pathways closely related to the occurrence and development of lung adenocarcinoma by analyzing gene expression microarrays [18]. In the present study, we used publicly available transcriptomic data and the mRNA expression-based stemness index (mRNAsi) as a quantitative reflector of cancer stemness to screen mRNAsirelated genes with prognostic potential. We then used the results to construct a risk score model for survival prediction in LUSC. Two independent cohorts were used to validate the prognostic performance of the risk score model. The molecular mechanisms underlying the survival subgroups were explored. The results of this analysis would contribute to the subtyping of survival groups and a more refined prognosis of LUSC.

Data source and retrieval
Gene expression data (FPKM value, Illumina HiSeq 2000 platform) from 501 tumor samples and 49 matched normal samples were downloaded from The Cancer Genome Atlas (TCGA) data repository (https:// gdc-portal. nci. nih. gov/). Among these, 494 tumor samples with corresponding clinical prognosis information were used as the training set (TCGA set).
We searched for the validation datasets in the NCBI GEO database using lung cancer and Homo sapiens as keywords. The inclusion criteria were as follows: histological information available, 150 or more samples, including 50 or more LSCC samples, and overall survival (OS) information of LSCC samples available. Consequently, the GSE30219 [19] (N = 307) and GSE37745 [20,21] (N = 196) datasets (GPL570 Affymetrix Human Genome U133 Plus 2.0 Array) met the criteria and were used as validation datasets in the current study, containing 58 and 66 LUSC samples, respectively, with corresponding clinical prognosis information.

Evaluation of mRNAsi scores and differentially expressed genes
Stem cell features of the tumor samples were evaluated using mRNAsi values, which were calculated using a oneclass logistic regression machine-learning algorithm in the gelnet package in R software (version 2.41-1); the protocol has been described in detail in our previous report [22].
We compared the mRNAsi scores of tumor and normal samples using the t-test in the R software program (version 3.6.1). Based on the median mRNAsi score, the tumor samples were categorized into low-and high-mRNAsi groups (values being below or above the median mRNAsi). Kaplan-Meier (KM) survival curves were plotted for each group and compared using the log-rank test.
We used the limma package [23] in R to screen for differentially expressed genes (DEGs), setting the threshold of significance at FDR < 0.05 when comparing tumor and normal samples and |log 2 FC|> 0.5 when comparing lowand high-mRNAsi samples. The genes that were common between the two lists of DEGs were then subjected to the gene ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using DAVID software (version 6.8, https:// david. ncifc rf. gov/). Statistical significance was set at P < 0.05.

Risk score model for survival prediction
Univariate Cox regression analysis was performed to identify survival-related common DEGs (log-rank p value < 0.05) using the survival package in R [24]. Of the survival-related genes, prognostic genes were identified by performing L1-penalized least absolute shrinkage and selection operator (LASSO) Cox regression analysis [25] using the penalized package. We used 1000-fold crossvalidation to determine the optimal lambda value, the penalty parameter, corresponding to the minimal meansquared error. Consequently, based on a linear combination of the LASSO Cox regression coefficients multiplied by the expression value of each optimal prognostic gene, the risk score was calculated for each sample in the TCGA set using the following formula: where β DEGs and Exp genes represent the regression coefficients and expression values, respectively.

Risk score PS) = β DEGs × Exp DEGs
Based on the median risk score, the tumor samples were then separated into high-and low-risk samples (having values above or below the median risk score).

Analysis of tumor-infiltrating immune cells and immune checkpoint genes
Differences in the fractions of tumor-infiltrating immune cells (TIICs) between high-and low-risk samples were analyzed using CIBERSORT [26] software. The expression levels of 14 immune checkpoint genes were also compared between the high-and low-risk samples.

Pathway enrichment analysis and protein prediction
Based on the gene expression data, we performed the KEGG pathway enrichment analysis using the Gene Set Enrichment Analysis described in the literature [27] (FDR < 0.05). We then accessed the Human Protein Atlas (HPV) database [28] (version 18, https:// www. prote inatl as. org/) to search for immunohistological images of the proteins encoded by the risk score genes in tumor and normal tissues.

Identification of mRNAsi-related DEGs
A flowchart depicting the overall design of the study is shown in Fig. 1. Using TCGA data, we found that mRNAsi was significantly higher in tumor tissues than in normal tissues (p < 0.001, Fig. 2a). Upon analyzing the median mRNAsi values, it was found that tumor samples could be categorized into high-and low-mRNAsi groups, with significantly different overall survival (OS) times (p = 0.00095, Fig. 2b). Additionally, age was significantly  Table 1).

Construction and validation of an eight-gene prognostic model
The results of the univariate Cox regression analysis indicate that a total of 77 mRNAsi-related DEGs were significantly associated with prognosis. Applying LASSO Cox regression analysis (1 mean squared error = 0.03871 (Additional file 1: Fig. S1)), we then selected a set of eight prognostic signature genes (PPP1R27, TLX2, ANKLE1, TIGD3, AMH, KCNK3, FLRT3, and PPBP). Based on the median expression level of each optimal signature gene, tumor samples were classified into high-and low-expression groups with significantly different OS times (p < 0.05, Fig. 4).
The expression data and the LASSO regression coefficients of the eight signature genes were used to calculate the risk score, as follows: The TCGA dataset was consequently dichotomized into high--and low-risk groups. In the ROC curve analysis, the 1-, 3-, and 5-year AUC values were 0.830, 0.749, and 0.749, respectively (Fig. 5a). The OS time was significantly longer in high-risk patients than in lowrisk patients (p = 0.00057, Fig. 6a). Furthermore, the eight-gene risk score model was applied to GSE37745 and GSE30219 datasets to validate the predictive performance of the model. The GSE37745 and GSE30219 datasets were divided by the eight-gene risk score into two risk groups with differential survival times (p = 0.029, 0.033, Fig. 6b, c; 1-year AUC = 0747, 0.783; In addition, we assessed the associations of mRNAsi levels with the risk score in the TCGA dataset. As shown in Fig. 7, the mRNAsi level was positively correlated with the risk score (PCC = 0.4402, p-value = 2.2E−16) and was significantly elevated in low-risk samples compared to high-risk samples (p-value = 8.1E−13).

Identification of independent prognostic factors and stratification analysis
Using the TCGA data, with clinical factors and risk score model status as variables, we performed univariate and multivariate Cox regression analyses to identify prognostic factors. In Table 3, we show that recurrence (HR = 2.047, 95%CI = 1.475-2.843, p-value = 1.880E−05) and risk score (HR = 1.571, 95%CI = 1.128-2.186, p-value = 7.530E−03) were independent prognostic factors.
As depicted in Fig. 8a, significantly different OS times were observed between patients with (N = 286) and without recurrence (N = 100, p < 0.0001). Stratification analysis was conducted according to recurrence. In patients without tumor recurrence, a worse prognosis was observed in the high-risk subgroup than in the low-risk subgroup (p = 0.0012, Fig. 8b). Regarding patients with tumor recurrence, the difference in OS time was insignificant between the two risk subgroups (p = 0.67, Fig. 8c).

The two risk groups had distinct immune characteristics and were significantly involved in DNA-repair-related pathways
There is evidence that cancer stemness is associated with immune checkpoint genes and the proportion of TIICs in the tumor microenvironment [22]. Therefore, we assayed the proportion of different types of TIICs. Compared to high-risk samples, low-risk samples had significantly increased percentages of naïve B cells (p = 0.006), CD8+ T cells (p = 0.044), and M1 macrophages (p = 0.022) and decreased percentages of resting memory CD4+ T cells (p = 0.022), monocytes (p = 0.01), and activated mast cells (p = 0.001, Fig. 9). We compared the expression of 18 immune checkpoint genes. Notably, CD47, HAVCR2, SIRPA, ICOS, TNFRSF9, BTLA, and TNFRSF4 were significantly downregulated in low-risk samples compared to high-risk samples (p < 0.05, Fig. 10). This result indicates that the two risk samples had different immune characteristics.
Moreover, we identified nine DNA-repair-related KEGG pathways that were significantly associated with the obtained risk subgroups (FDR < 0.05), including base excision repair, RNA degradation, RNA polymerase, and spliceosome (Table 4). Using data from the HPA database, we found immunohistochemical images of five prognostic signature genes in normal and tumor tissues (Fig. 11), including three upregulated genes (ANKLE1, PPP1R27, and AMH) and two downregulated genes (FLRT3 and PPBP). The immunohistochemical results were consistent with our differential expression analysis (Fig. 11, Additional file 2: Table S1).

Discussion
Functionally defined by their high tumorigenic potency and self-renewal properties, CSCs are a critical driving force for cancer metastasis, recurrence, and chemoresistance and have been increasingly acknowledged as potential therapeutic targets [29]. In the present study, we focused on the identification of CSC-related prognostic genes for survival prediction in LUSC. By exploiting gene expression data from TCGA, we identified 404 mRNAsi-related DEGs in tumors-77 of which were significantly associated with survival-and constructed a risk score model using eight prognostic genes obtained using LASSO Cox regression analysis. The eight-gene risk score partitioned the TCGA dataset into two risk groups with significantly different OS times (p = 0.00057, 1-, 3-, and 5-year AUC = 0.830, 0.749 and 0.749, respectively). The value of mRNAsi was positively correlated with the risk score. The capability of the eight-gene risk score to stratify survival into subgroups was successfully validated in the two validation datasets, as evidenced by the significant log-rank p-values and high AUC values (0.7-0.8). Furthermore, the eight-gene risk score was shown to be an independent prognostic factor, regardless of recurrence rate. These results extend our knowledge of the maintenance and promotion of the malignant characteristics of CSCs and may contribute to a more accurate prognosis (survival prediction) as well as targeted therapies for LUSC. The eight-gene prognostic signature identified in our study was composed of PPP1R27, TLX2, ANKLE1, TIGD3, AMH, KCNK3, FLRT3, and PPBP. Anti-Mullerian hormone (AMH), a member of the transforming growth factor/bone morphogenetic protein superfamily, participates in regulating epithelial-mesenchymal transition (EMT), epithelial plasticity, and chemoresistance in lung cancer [30]. Moreover, AMH is an immune-related prognostic gene in LUSC and has been used to construct a prognostic model in LUSC [31]. Zhuang et al. also found that AMH is a LUSC-related immune gene and is not associated with distant metastasis [32]. TWIK-related acid-sensitive potassium channel 1 (TASK1), encoded by KCNK3, is associated with pulmonary circulation and controls pulmonary arterial tone, which may contribute to poor prognosis in lung cancer patients [33]. KCNK3 has been shown to influence apoptosis and proliferation in NSCLCs, and KCNK3 knockdown increases apoptosis in tumor cells [34]. Therefore, KCNK3 may play a role in cell motility, activation and proliferation [35]. FLRT3 is a transmembrane protein belonging to the family of axon guidance-related factors. FLRT3, which is found in many tissues and is involved in cell adhesion and adipocytokine signalling pathways [35], has been implicated in the progression and prognosis of LUSC and could serve as a prognostic biomarker [36]. Pro-platelet basic protein (PPBP) and chemokine ligand 7 (CXCL7) are platelet activation markers that act as inducers of macrophage chemotaxis and mediators of neutrophil accumulation [37]. PPBP is a survival-related hub gene in lung adenocarcinoma [18] and non-smoker females with lung cancer [38]. Furthermore, studies have shown that PPBP is significantly increased in lung cancer tissue and blood samples, making it a novel diagnostic biomarker for lung carcinoma [39,40]. T-cell leukemia homeobox 2 (TLX2) has prognostic value in uterine sarcoma [41]. ANKLE1 (ankyrin repeat and LEM domain) is involved in DNA damage response and DNA repair and is associated with breast cancer development [42,43]. In addition, ANKLE1 has been repeatedly shown to be involved in DNA repair pathways in preclinical and in vitro screens, including endonuclease activity, proliferation, and drug response in CRISPR screens of cancer cell lines [44][45][46]. A study of NSCLC showed that ANKLE1 RNAi in combination with paclitaxel increased the efficacy of drug response [47]. Nonetheless, there is little information regarding the biological functions of PPP1R27 and TIGD3 in cancer. In this study, the capability of the eight-gene risk score to stratify survival time was successfully validated using two validation datasets, suggesting that these eight genes may be prognostic biomarkers for LUSC.
The immune environment is predictive of the prognosis of NSCLCs [48]. TTIICs have been implicated extensively in the initiation and progression of LUSC [49]. In the current study, increased numbers of CD8+ T cells and M1 macrophages and significant downregulation of immune checkpoint genes were observed in low-risk patients compared to high-risk patients. These results showed that the high -and low-risk subgroups possessed distinct immune microenvironment characteristics. CD8+ T cells mediate an anti-tumor immune response [50], and M1 macrophages exert pro-inflammatory and anti-tumor actions [51]. Immune checkpoints are critical for immune suppression and evasion in cancers, and immune checkpoint inhibitors represent an efficient therapeutic approach against a wide spectrum of malignancies [52]. Our results suggest that low-risk patients have a stronger antitumor immune function, which protects against LUSC and achieves a significant survival benefit. Through interactions with the tumor immune microenvironment, CSCs facilitate immune evasion and suppress the immune system to promote tumor progression [53]. Taken together, these results indicate that the crosstalk between CSCs and the immune microenvironment may affect the prognosis of LUSC patients, and further studies are needed to validate these results. This study has some limitations. Firstly, because all data were obtained from the TCGA and GEO databases, selection bias could not be ruled out. Secondly, 501 LUSC tumor samples and 49 normal control samples were downloaded from the TCGA database, considering that such unequal sample distribution (controls being approximately 10% of tumor samples), may amplify the detection of differences. Thus, further validation is required to support the discovery of this research.

Conclusion
In this study, we constructed and validated a risk score model based on the expression of eight CSC-related DEGs that could effectively predict LUSC outcomes. These eight CSC-related genes may be prognostic biomarkers and potential therapeutic targets for LUSC. Our study sheds light on the prognostic value of cancer stemness-related genes and their underlying mechanisms and may facilitate personalized counselling and treatment of LUSC. Further research is required to confirm and extend these findings.